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INTRODUCTION 


The  long  term  objective  of  the  present  effort  is  the  development  of  solution 
techniques  for  direct  numerical  simulation  of  unsteady  3-D  incompressible  turbulent 
flows.  The  kinetic  aspects  of  this  problem  are  governed  by  a  set  of  parabolic  partial 
differential  equations,  which  may  be  efficiently  integrated  by  a  variety  of  time  marching 
schemes.  The  kinematic  aspects  of  this  flow  such  as  the  relationship  between  velocity 
and  vorticity,  and  the  relationship  between  velocity  and  pressure  are  governed  by  elliptic 
partial  differential  equations,  which  can  be  solved  at  any  instance  in  time,  only  by 
iterative  techniques.  Direct  and/or  large  eddy  simulation  of  turbulent  flows  over 
submarine  configurations,  turbomachinery,  pumps,  ducts  and  other  configurations  of 
interest  to  the  U.  S.  Navy  require  efficient  solution  methods  for  solving  the  governing 
equations. 

The  near  term  objective  of  the  present  research  is  to  investigate  and  develop 
efficient  time  marching  schemes  for  integrating  the  governing  equations,  and  to  evaluate 
the  stability  and  accuracy  of  the  schemes  developed  by  studying  a  class  of  2-D  and  3-D 
unsteady  external  flows  for  which  good  quality  experimental  and  analytical  results  are 
available. 


WORK  DONE  DURING  THE  REPORTING  PERIOD 

During  the  reporting  period,  extension  of  the  2-D  unsteady,  incompressible 
viscous  flow  methodology  was  to  three-dimensions  was  completed.  A  very  general  3-D, 
incompressible  flow  solver  capable  of  handling  arbitrary  curvilinear  grids  has  been 
developed.  The  grid  may  move  or  deform  with  time,  as  will  be  the  case,  for  example,  for 
incompressible  viscous  flow  past  a  spinning  propeller.  The  scheme  is  third  order  accurate 
in  space,  and  first  or  second  order  accurate  in  time. 

This  solver  has  been  applied  to  the  following  cases: 

a)  Incompressible  viscous  flow  past  an  ellipsoid  at  an  angle  of  attack.  This 
geometry  was  chosen  because  of  the  availability  of  existing  experimental  data. 


b)  Incompressible  viscous  flow  through  a  90  degree  bend,  of  rectangular  cress 
section. 


c)  A  multi-grid  scheme  has  also  been  implemented  for  acceleration  of  the 
iterative  process  at  every  time  step.  Preliminary  results  on  a  two  grid  sequence  are 
encouraging. 

The  above  results  are  discussed  in  detail,  along  with  the  complete  mathematical 
and  numerical  formulation  in  the  following  two  documents: 

1.  W.  G.  Pork,  "Numerical  Solution  of  3-D  Unsteady  Incompressible 
Viscous  Flows,"  Ph.  D.  Thesis  Proposal,  to  be  presented  to  the  Thesis  Committee  in 
August  1992. 

2.  W.  G.  park  and  L.  N.  Sankar,  "  A  Technique  for  the  Prediction  of 
Unsteady  Incompressible  Viscous  Flows,"  Abstract  submitted  to  the  AIAA  Aerospace 
Sciences  Meeting  in  Reno,  Nevada  1992. 

Copies  of  the  above  two  documents  are  enclosed. 

ANTICIPATED  RESULTS  FOR  THE  NEXT  REPORTING  PERIOD 

By  the  conclusion  of  the  next  reporting  period  (November  30,  1991),  we  plan  to 
have  the  following  case  completed: 

a)  Viscous  flow  over  a  highly  twisted  tapered  spinning  propeller  in  forward 
motion.  An  SR-7  propeller  (developed  for  aircraft  applications)  is  being  used  for  the 
code  validation.  But  the  flow  solver  can  handle  marine  propeller  configurations  as  well. 

TECHNOLOGY  TRANSFER 

Dr.  Wei  Tang  of  the  Naval  Research  center  at  Annapolis,  Maryland  (Phone:  301- 
267-2730)  has  acquired  a  version  of  our  2-D  unsteady  incompressible  flow  solver,  and 
plans  to  extend  it  to  flow  through  multi-stage  turbine  and  compressor  configurations.  We 
hope  to  assist  her  in  the  validation  of  the  flow  solver,  as  needed. 
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CHAPTER  I 


INTRODUCTION 


The  accurate  computation  of  three-dimensional  unsteady  incompressible  flow 
problem  is  one  of  great  interest  to  researchers  working  in  fields  of  aerodynamics, 
hydrodynamics  and  biofluid  mechanics.  The  flow  over  complex  submarine  shapes,  flow 
past  underwater  propeller,  flow  within  turbomachinery,  and  flow  in  blood  vessels  with 
compliant  walls  are  examples  of  such  flows.  Accurate  and  efficient  computation  of  such 
flows  at  high  Reynolds  numbers  is  presently  not  possible  due  to  the  mixed  (elliptic- 
parabolic)  nature  of  the  governing  equations.  Indeed,  methods  for  three-dimensional 
incompressible  flows  lag  behind  three-dimensional  compressible  flows  by  several  years. 
Until  accurate  and  efficient  methods  for  three-dimensional  incompressible,  unsteady  flows 
become  available,  it  will  not  be  possible  to  attempt  challenging  problems  such  as  the  first 
principles  based  on  direct  simulation  or  large  eddy  simulation  of  turbulent  flow  s  over 
complex  geometries.  The  lack  of  such  tools  is  one  of  the  principal  reasons  that  the  first 
principles  based  prediction  of  turbulent  flows  past  and  through  complex  configurations  has 
not  been  extensively  attempted  to  date. 

As  Gresho  and  Sani  (ref.l)  pointed  out,  the  pressure  is  a  somewhat  mysterious 
quantity  in  incompressible  flows.  It  is  not  a  thermodynamic  variable  since  there  is  no 
'equation  of  state'  for  an  incompressible  fluid.  It  is  in  one  sense  a  mathematical  artefact  -  a 
Lagrange  multiplier  that  constrains  the  velocity  field  to  remain  divergence-free  ;  i.e. 
incompressible  -  yet  its  gradient  is  a  relevant  physical  quantity  ;  a  force  per  unit  volume.  It 
propagates  at  infinite  speed  in  order  to  keep  the  flow  always  and  everywhere 
incompressible  ;  i.e.  it  is  always  in  equilibrium  with  a  time-varying  divergence-free  velocity 
field. 


One  might  have  the  idea  that  the  compressible  Navier-Stokes  equation  solvers  can 
compute  incompressible  flows  using  compressible  flow  methods,  and  setting  the  Mach 
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number  to  be  very  low.  But  this  idea  becomes  impractical  at  very  low  Mach  numbers 
because  the  compressible  Navier-Stokes  equation  solvers  have  a  singular  behavior  as  the 
Mach  number  approaches  zero.  This  leads  to  an  ill-conditioned  stiff  system  of  equations 
and  consequently  very  slow  convergence,  or  even  divergence  of  the  solution  with  time. 
.This  stiffness  can  be  explained  as  a  time  step  limitation  (ref.2).  We  note  that  all  explicit 
-methods  for  solving  the  compressible  Navier-Stokes  equations  are  limited  to  a  time  step 
*  which  is  less  than  that  given  by  the  CFL  condition.  For  example,  in  two-dimensions  : 

At< _ 1 _ 

(|u| /  Ax)  +  (jv|  /  Ay)  +  a  ^(1  /  Ax)2  +  (1  /  Ay)2 j  * 

where  a  is  the  speed  of  sound.  From  this  condition,  we  observe  that  At  approaches  zero 
as  the  speed  of  sound  approaches  infinity.  As  a  result,  an  "infinite"  amount  of  computer 
time  would  be  required  to  compute  a  truly  incompressible  flow  in  this  manner.  Implicit 
methods  will  permit  a  larger  At ,  but  the  maximum  value  is  normally  less  than  100  times 
that  given  by  Eq.(l.l)  because  of  truncation  errors,  approximate  factorization  errors,  and 
so  on.  Thus,  even  if  an  implicit  scheme  is  used,  it  is  not  practical  to  compute  a  truly 
incompressible  Navier-Stokes  solution  using  compressible  flow  methods. 

The  significant  difficulty  in  solving  incompressible  Navier-Stokes  equations  is  that 
the  governing  equations  are  a  mixed  elliptic-parabolic  type  of  partial  differential  equations 
The  continuity  equation  does  not  have  a  time  derivative  term  and  is  given  in  the  form  of  a 
divergence-free  constraint.  This  is  another  major  difference  between  the  incompressible  and 
compressible  Navier-Stokes  equations.  The  absence  of  a  time  derivative  term  in  the 
continuity  equation  prohibits  time  integration  of  continuity  equation  by  a  time  marching 
scheme.  The  compressible  Navier-Stokes  equations,  on  the  other  hand,are  efficiently 
integrated  by  time  marching  schemes  because  they  are  a  set  of  parabolic  partial  differential 
equations. 

One  of  the  commonly  used  approaches  for  solving  two-dimensional  incompressible 
flow  is  the  vorticity-velocity  or  vorticity-stream  function  formulation  (ref.  3,4,5).  This  is 
very  efficient  for  two-dimensional  problems,  but  this  approach  can  not  be  extended 
straightforwardly  to  three  dimensions.  Consequently,  the  incompressible  Navier-Stokes 
equations  for  three-dimensional  problem  are  normally  solved  in  their  primitive  variable 
form  (p,u,v,w).  Most  methods  using  primitive  variables  may  be  classified  into  three 
groups.  The  first  approach  is  the  pressure  Poisson  method  or  Marker-and-Cell  (MAC) 


method  which  was  first  introduced  by  Harlow  and  Welch  (ref.6).  In  the  pressure  Poisson 
method,  the  velocity  field  is  advanced  in  time  by  solving  the  momentum  equations  with  a 
stable  explicit  or  implicit  time  marching  scheme.  Then  the  pressure  field  is  evaluated  at  each 
time  step  by  solving  a  Poisson  equation  for  pressure  directly  (ref.7)  or  iteratively 
.(ref. 8, 9, 10).  The  continuity  equation  is  thus  satisfied  when  the  pressure  field  is  computed 
implicitly.  This  Poisson  equation  for  pressure  is  obtained  by  taking  the  divergence  of  the 
‘unsteady  momentum  equations.  The  main  idea  of  the  MAC  method  (ref.11,12),  an 
alternative  to  solving  a  pressure  Poisson  equation,  is  that  the  pressure  field  is  updated  at 
each  time  step  by  adjusting  the  pressure  by  an  amount  proportional  to  the  negative  of  the 
velocity  divergence : 

Pjj-Pfj1  — P(V-V)k-1  (1.2) 

Here  the  superscripts  ’k’  and  'k-l*  denote  the  iteration  level,  and  (3  is  a  relaxation  factor. 
Usually,  a  staggered  grid  system  (ref.6)  is  used  for  the  MAC  method,  because  such  a  grid 
does  not  require  the  specification  of  pressure  on  the  boundaries  and  does  not  produce 
unphysical  oscillations  in  the  pressure  and  velocity  fields  due  to  the  central  differencing  of 
the  pressure  gradient  term.  The  second  approach  is  a  projection  method  (or,  fractional  step 
method  )  which  was  first  introduced  by  Chorin  (ref.  13).  At  the  first  step,  an  intermediate 
velocity  is  computed  from  the  momentum  equation  without  the  pressure  gradient  term. 
Then  a  pressure  field  is  computed  which  will  make  the  velocity  field  obtained  from  the  first 
fractional  step  divergence  free.  Finally,  a  second  fractional  step  is  performed  using  the 
pressure  field  just  computed.  The  third  group  is  the  pseudocompressibility  method 
(ref.14,15)  which  was  also  first  introduced  by  Chorin  (ref.  1 6)  primarily  for  obtaining 
steady  state  solutions.  In  this  method,  an  artificial  pressure  derivative  with  respect  to  time  is 
appended  to  the  continuity  equation.  The  entire  system  of  equation  is  solved  by  a  time 
marching  scheme  developed  for  compressible  flows,  such  as  the  approximate  factorization 
scheme  (ref.  17).  If  only  a  steady  state  is  of  interest,  then  the  added  pressure  derivative 
drops  out  in  the  steady  state,  and  physically  correct  solutions  are  achieved.  If  the  aim  is  to 
achieve  time-accurate  calculations,  either  the  artificial  pressure  derivative  should  be  kept 
very  small  (which  makes  the  equations  extremely  stiff,  and  forces  very  small  time  steps)  or 
an  inner  iterative  loop  within  each  time  step  should  be  used  (ref.  18,19).  A  concept  similar 
to  the  pseudocompressibility  method,  known  as  the  penalty  function  method  (ref.20)  is 
widely  used  in  the  finite-element  based  incompressible  flow  solvers,  which  solves  for  p  to 
satisfy  : 
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Xp  +  VV=0  (1.3) 

In  this  method,  the  pressure  gradient  term  of  momentum  equation  is  eliminated  by 
substituting  Eq.(1.3)  into  the  momentum  equation,  and  then  solving  the  momentum 
equations  with  X  — >  0. 

The  methods  for  solving  incompressible  viscous  flow  discussed  above  have  several 
drawbacks : 

a)  Most  of  them  are  only  second  order  accurate  in  space,  and  first  or  second  order  accurate 
in  time.  Before  these  schemes  can  be  applied  to  phenomena  such  as  direct  numerical 
simulation  of  turbulence,  it  will  be  necessary  to  raise  the  spatial  and  temporal  accuracy  to 
fourth  or  higher  order. 

b)  The  iterative  convergence  of  the  pressure  Poisson  solvers  deteriorates  at  high  Reynolds 
numbers. 

c)  In  some  instance  (e.g.  in  the  pseudocompressibility  method),  a  trade  off  exists  between 
temporal  accuracy  and  convergence  speed. 

d)  These  methods  do  not  take  advantage  of  the  vast  progress  that  has  been  achieved  in  the 
solution  of  steady,  viscous  flows.  For  example,  with  rare  exceptions,  multigrid 
acceleration  of  Poisson  solvers  has  not  been  attempted.  Acceleration  of  the  iterative  solution 
of  the  pressure  field  to  convergence  using- spatially  varying  time  steps  and  grid  sequencing 
have  also  not  been  extensively  used. 

e)  There  has  been  a  growing  interest  in  the  use  of  massively  parallel  computer  architectures 
such  as  the  Connection  Machine  to  solve  unsteady  viscous  flows.  Many  of  the 
compressible  flow  algorithms  have  already  been  adapted  for  use  on  these  machines.  There 
is  a  need  to  develop  new  procedures  and  modify  existing  algorithms  for  incompressible 
flows,  on  parallel  machines. 

The  objective  of  this  study  is  to  develop  an  efficient  and  accurate  solution 
technique  for  the  analysis  of  two-  and  three-dimensional,  unsteady,  incompressible, 
viscous  flows.  The  key  features  of  the  present  scheme  are  listed  below  : 

a)  The  primitive  variables  (p,u,v,w)  are  the  primary  unknowns  in  the  present  formulation. 

b)  The  equations  and  the  solution  procedures  are  cast  into  a  curvilinear,  time-deforming 
coordinate  system  to  handle  complex  internal  and  external  flows. 

c)  An  iterative  time-marching  scheme  is  used. 

d)  The  present  scheme  is  semi-implicit  at  each  iteration  and  is  suitable  for  efficient 
execution  on  the  current  generation  of  vector  or  massively  parallel  computer  architectures. 


c)  The  solution  procedure  works  for  a  wide  range  of  Reynolds  numbers,  with  no 
appreciable  loss  in  solution  efficiency. 

f)  The  present  scheme  is  first  order  accurate  in  time  and  second  order  accurate  in  space,  but 
higher  order  accuracy  in  space  and  time  is  easily  achievable. 

Only  laminar  flow  is  considered  in  the  results  to  be  discussed  because  the  goal  of 
-this  study  is  to  develop  an  efficient  and  accurate  incompressible  Navier-Stokes  solver.  This 
‘method  is  however  capable  of  handling  turbulent  flows  provided  a  suitable  turbulence 
model  is  used,  and  there  are  no  inherent  limitations  in  this  method  that  will  restrict  it  to 
laminar  flows. 
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CHAPTER  II 

MATHEMATICAL  FORMULATION 


In  this  chapter,  the  governing  equation  for  three-dimensional,  unsteady, 
incompressible,  viscous  flow  are  presented  in  terms  of  the  primitive  variables  (p,u,v,w)  in 
both  the  Cartesian  coordinate  system  and  a  curvilinear  non-orthogonal,  time  deforming 
coordinate  system. 

2.1  Governing  Equations  in  the  Physical  Domain 

The  motion  of  an  incompressible  viscous  flow  is  governed  by  the  conservation  of 
mass  and  momentum,  so  called  the  continuity  equation  and  the  Navier-Stokes  equation. 
Three-dimensional  unsteady,  incompressible,  laminar,  Navier-Stokes  equations  in  an 
inertial  Cartesian  coordinate  system  may  be  written  in  a  non-dimensional  form  as  follows  : 

|  +  |r(E-E>)+|r(F-F.)  +  ^(G-O,)  =  0  (2.1, 

where 

q 


The  stress  terms  are  given  by 
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In  Eq.(2.2)  and  Eq.(2.3),  u,  v  and  w  are  the  normalized  Cartesian  components  of  velocity, 
p  is  the  normalized  pressure,  and  Re  is  the  Reynolds  number  defined  as  : 


(2.4) 


where  p,  L  and  p  are  fluid  density,  freestream  velocity,  reference  length  and 
coefficient  of  viscosity  (dynamic  viscosity),  respectively. 

The  governing  equation  (2.1)  is  a  mixed  set  of  ell’~ric -parabolic  partial  differential 
equations.  As  mentioned  before,  the  absence  of  a  time  der'  ative  in  continuity  equation  and 
the  absence  of  an  explicit  relationship  between  pressure  and  divergence-free  condition  on 
the  velocity  prohibit  time  integration  in  a  straightforward  manner  by  a  stable  time  marching 
scheme.  In  this  study,  the  continuity  equation  is  modified  to  directly  link  the  iterative 
changes  in  pressure  to  changes  in  velocity,  as  done  in  the  Marker-and-Cell  method. 


2.2  Governing  Equations  in  the  Computational  Domain 

If  the  above  equations  are  directly  used  on  a  Cartesian  system  to  flow  past  complex 
geometries,  the  imposition  of  boundary  conditions  will  require  a  complicated  interpolation 
of  the  data  on  local  grid  lines,  since  the  computational  boundaries  of  complex  geometries 
do  not  coincide  with  coordinate  lines.This  leads  to  a  local  loss  of  accuracy  in  the  computed 
solution  and  leads  to  a  complex  program.  To  avoid  these  difficulties,  a  transformation  from 
the  physical  domain  (Cartesian  coordinates! t,x,y,z))  to  computational  domain  (generalized 
curvilinear  coordinates(t,^,r|,^))  is  used.  After  transformation  from  the  physical  domain  to 
the  computational  domain,  the  governing  equations  can  be  written  as  : 
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1 
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(VC-V^+(VC.Vii)v,  +  (VC.VC)v5 
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(2.5) 


(2.6) 


(2.7) 


with  the  contravariant  velocities  U,  V  and  W  : 


U  =  +  U^X  +  V^y  +  W^Z 

V  =  Ht  +  urix  +  vTiy  +  wrjz  (2.8) 

W  =  C  +  uCx  +  vCy  +  WC 


Here  J  is  the  Jacobian  of  transformation 
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j  _  gfenj)  _  i 

9(x,y,z) 

y%  y;  (2.9) 

Z5  zn  Z? 

-The  quantities  ,  qt  and  316  presented  if  the  grid  is  in  motion  (as  in  the  case  of  flow 
past  an  oscillating  airfoil  or  a  spinning  propeller).  These  quantities  are  given  in  terms  of  the 
velocity  of  the  grid  (xt,  yt,  zT)  with  reference  to  a  stationary  observers  : 

Ct  =  —  Xx  Cx  —  y-t  ^y  —  zx  Cz 
"Ht  =  —  Xx  ^lx  —  yx  ^ly  —  zx  *lz 
Cl  =  -  XX  Cx  -  yx  Cy  -  ZX  Cz 


(2.10) 


CHAPTER  III 


1  0 


NUMERICAL  FORMULATION 


The  numerical  procedure  for  solving  the  governing  equation  is  an  iterative  time 
marching  scheme  which  attempts  to  solve  the  discretized  form  of  equations  to  a  user- 
specified  accuracy  at  any  time  step.  Details  of  the  iterative  process  are  given  in  this  chapter. 

3.1  Grid  Generation 

The  present  method  is  a  finite  difference  scheme  which  solves  the  discretized  form 
of  the  partial  differential  equations  at  a  set  of  discrete  points  in  the  flow  field.  Therefore,  a 
set  of  grid  points  within  the  domain,  including  its  boundaries,  must  be  specified  before 
solving  the  governing  equations.  Such  a  body-fitted  grid  system  may  be  generated  by 
conformal  mapping,  by  algebraic  method,  or  by  partial  differential  equation  techniques.  In 
this  study,  body-fitted  C-grid  (Fig.l)  and  H-O  grid  system  (Fig.5)  are  generated  by  an 
algebraic  method  for  two-dimensional  flow  around  NACA  0012  airfoil  and  three- 
dimensional  flow  around  the  ellipsoid  of  revolution,  respectively.  For  the  three- 
dimensional  curved  duct  problem,  a  sheared/rotated  Cartesian  grid  is  used. 

3.2  Grid  Motion 

In  unsteady  state  computations,  it  is  convenient  to  use  a  moving  grid  to  account  for 
the  body  motion.  The  grid  is  attached  to  the  body  and  it  rotates  or  translate  with  the  body. 
The  grid  coordinates  can  be  advanced  explicitly  by  a  first  order  time  marching  scheme  : 


xn+1  =  x"  +  x"  At 
yn+1  =  yn  +  y"  At 

zn+l  =  zn  +  zn  fa 


(3.1) 


However,  if  only  a  pure  rotational  motion  is  considered  (say  in  a  two-dimensional  flow 
problem),  new  coordinates  of  grid  at  any  instance  in  time  can  be  simply  obtained  by  using 
the  following  relations : 


x 
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’cos0 

-sin0 

y 

sin0 

cos  0_ 

z' 

(3.2)  . 


where  (x,  z)  is  the  instantaneous  x,  z  values  of  the  node  and  (x',z')  is  the  x,  z  values  of  the 
node  prior  to  rotation,  and  0  is  the  clockwise  rotation  angle.  In  such  a  case  xT  and  zt  may 
be  found  by  analytical  differentiation  of  (3.2)  with  respect  to  time  or  from  (3.1). 

3.3  Iterative  Time  Marching  Procedure 

The  goal  of  the  present  procedure  is  to  advance  the  flow  properties  (p,u,v,w)  from 
a  known  time  level  'n'  to  the  next  time  level  ‘n+l’.  First  of  all,  let  us  consider  the 
momentum  equation.  Since  the  momentum  equation  is  a  parabolic  type  of  partial  differential 
equation,  it  can  be  solved  using  a  time  marching  scheme  as  follows  : 


_  5-)  +  8sE"*“  +  8 ,r'm  +  8(C  ™  = 

5tEn;"+5,F';’+stcr 

where  q  is  §  of  Eq.(2. 6)  excluding  the  first  row  element,  i.e., 

u 
v 
w 


(3.3) 


(3.4) 


Similarly,  E,  F,  G,  Ev,  Fvand  Gv  can  be  also  defined.  For  example. 


E  = 


uU  +  p4x 
VU  +  P^y 


wU  +  p^J 


(3.5) 


The  above  discretization  of  Eq(3.3)  is  first  order  accurate  in  time  if ’m'  is  zero  or  one,  and 
second  order  accurate  if 'm'  is  set  to  1/2.  The  operators,  5^,5^  and  8^  represent  second 

order  accurate  or  higher  order  accurate  spatial  differences.  The  higher  order  spatial  accuracy 
may  be  achieved  on  uniform  grids  using  Pade  approximations  to  the  derivatives;  on  highly 
-stretched  grids,  higher  order  accuracy  may  be  achieved  using  a  Lagrangean  fit  to  the  flow 
'variables.  In  high  Reynolds  number  flows,  the  Lagrangean  fit  need  not  be  equally  weighted 
about  the  node,  but  may  be  biased  in  the  direction  of  flow.  For  example,  when  the  flow  is 
from  left  to  right,  if  the  Lagrangean  interpolation  of  flow  variables  is  done  using  nodes 
only  to  the  left  of,  and  including,  the  current  node,  then  an  upwind  formulation  results. 

If  the  Newton  iteration  method  is  applied  to  solve  this  unsteady  flow  problem, 
Eq.(3.3)  is  rewritten  as  follows  : 


JL(qn+1- k+1 

At 


-  qn)  +  5^En+m’ k+1  +  8T1Fn+m- k+1  +  6;G 


/-'•n+m,  k+1  _ 


6rn+m,  k+1  ,  c  k+l  ,  K  7Fn+l,  k+1 

j;Ev  +  S^F,,  +  8^GV 


(3.6) 


Following  a  local  linearization  of  E,  F,  G,  Ev,  F  and  Gv  about  the  'n+m'  time  level  and 
at  the  'k'  iteration  level,  one  may  have 


I  +  —  A  +  —  B  +  —  cl  Aq  =  0)Rn+m’k 

{  *  an  ac  J 


where  to  is  a  relaxation  factor  and  A,  B  and  C  are  the  Jacobian  matrices  of  the  flux  vectors 
E  -  Ev,  F  -  Fv  and  G  -  Gv ,  respectively: 


9(E-EV) 

9q 


B  = 


a(F-Fv) 

9q 


C  = 


8(G  -  Gv) 

dq 


(3.8) 


and  Rn+m>  k  is  the  residual  vector,  defined  as : 


-n+l,  k  -n 
tn+m,  k  _  0 _ Q _ 


At 


-  (5,En+m- k  +  511Fn+mi  k  +  5;-Gn+rn'  k ) 


ihE 


Fn+m,  k  ,  c  r-n+m,  k  ,  j*  n  n+m,  k 
v  +6nFv  +o^Gv 


) 


(3.9) 


Note  that  when  Rn+m’  k  goes  to  zero,  the  momentum  equations  in  their  discretized  form  are 
exactly  satisfied,  and  the  solution  is  independent  of  co,  and  any  approximations  made  in  the 
construction  of  A,  B  and  C. 

Next,  let's  consider  the  continuity  equation.  As  mentioned  in  Chapter  I,  in  order  to 
.solve  incompressible  viscous  flow  problems  efficiently,  we  need  a  relationship  coupling 
-changes  in  the  velocity  field  with  changes  in  the  pressure  field  while  satisfying  the 
‘divergence -free  constraint.  In  the  present  study,  the  Marker-and-Cell  (MAC)  approach  is 
used  to  link  the  iterative  changes  between  them,  and  can  be  written  : 

Ap  =  -p  (V  ■  V)n+1‘ k  (3.10) 


,  A _ n+1,  k+1  _n+l,  k 

where  Ap  =  p  -  p 

and  p  is  a  relaxation  factor,that  may  even  vary  from  node  to  node  using  local  time  concept. 
Again,  when  Ap  goes  to  zero,  the  continuity  equation  is  exactly  satisfied  at  each  time  step, 
even  in  unsteady  flows. 

In  curvilinear  coordinate  system,  Eq.(3.10)  can  be  rewritten  as  : 
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(3.11) 


The  contravariant  velocities,  U,  V  and  W  are  already  defined  in  Eq.(2.8). 

Eq.(3.10)  states  that  if  a  cell  is  accumulating  mass,  then  the  pressure  value  at  next 
iteration  is  increased  to  repel  fluid  away  from  the  cell.  If  a  cell  is  losing  mass,  then  the 
pressure  value  is  lowered  to  draw  fluid.  Thus  the  pressure  field  is  iteratively  updated  along 
with  the  velocity  field  until  the  conservation  of  mass  is  satisfied. 

Combining  the  momentum  equation,  Eq.(3.7)  and  the  continuity  equation, 
Eq.(3.1 1) ,  and  applying  the  numerical  discretization  in  time  and  space  at  all  nodes  in  the 
flow  field,  a  system  of  simultaneous  equation  results  for  the  quantity  A(j  equal  to 


D  U  V  W  A 

Ay,  Ay,  Ay,  AyJ.  This  system  may  be  formally  written  as  : 

[M]{A§}  =  (R} 


(3.12) 


Here,  since  the  right  hand  side  is  the  discretized  form  of  the  unsteady  governing 
equations,  as  long  as  is  driven  to  zero,  the  discretized  form  of  unsteady  Navier- 

Stokes  equations  are  exactly  satisfied  at  physical  time  level  'n+1'. 

Although  the  matrix  [M]  is  a  sparse,  banded  matrix,  direct  inversion  of  this  matrix 
-requires  a  huge  number  of  arithmetic  operations.  A  common  strategy  in  iterative  solutions 
"of  elliptic  equations  is  to  approximate  the  matrix  [M]  by  another,  easily  inverted  matrix 
*  [N] .  The  closer  the  matrix  [N]  is  to  [M] ,  the  faster  the  iterative  convergence  of  the 
solution  at  any  time  step.  In  this  study,  matrix  [N]  contains  only  the  diagonal 
contributions  of  matrix  [M] ,  and  Eq.(3.12)  becomes  an’ explicit  form  which  is  easier  to  be 
tailored  for  efficient  execution  on  the  current  generation  of  vector  or  massively  parallel 
computer  architectures  than  an  implicit  form.  This  simplicity  comes  at  the  expense  of  the 
iterative  speed.  Acceleration  of  the  iterative  process  above  is  a  major  contribution  of  this 
work  to  the  state  of  the  art. 

The  spatial  derivatives  of  convective  flux  terms  are  differenced  by  using  third  order 
accurate  upwind  QUICK  (Quadratic  Upstream  Interpolation  for  Convective  Kinematics, 
ref.21)  scheme  to  reduce  unphysical  oscillations  or  false  diffusion  for  high  Reynolds 
number  flows,  and  the  spatial  derivatives  of  viscous  terms  are  differenced  using  half-point 
central  differencing.  The  spatial  derivatives  of  continuity  equation  is  differenced  with 
central  differencing  and  a  fourth  order  artificial  damping  term  is  added  to  the  continuity 
equation  to  stabilize  the  present  procedure.  The  QUICK  scheme  is  constructed  that,  instead 
of  such  a  linear  interpolation  for  the  convective  terms  as  used  in  standard  one-sided 
differencing  schemes,  a  three-point  upstream  weighted  quadratic  interpolation  is  used.  For 
example,  let’s  consider  the  convective  term  in  £ -direction  which  may  be  approximated  as 
follows : 


where 


uin  =  j_  ruin  ruin 

dt,{  J  J~A$  [l  J  Ll  l  J  Ji-i _ 
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(3.13) 


(3.14) 


The  curvature  terms  (CURV)  depend  on  the  direction  of  the  contravariant  velocity  U  : 
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(3:15) 


(3.16) 


(a) 


Fig.  3.1.  Quadratic  upstream  interpolation 

(a)  For  U  >  0 

(b)  For  U  <  0 


3.4  Initial  and  Boundary  Conditions 

The  governing  equation  (2.1)  and  (2.5)  is  a  mixed  elliptic-parabolic  type  of  partial 
differential  equation,  and  requires  initial  conditions  to  start  the  calculation  as  well  as 


boundary  condition  at  every  time  step.  The  parabolic  nature  of  the  flow  ensures  that  the 
flows  will  be  independent  of  initial  conditions,  after  large  number  of  time  step. 

In  the  present  work,  the  quantities  Ap,  Au,  Av  and  Aw  are  set  to  zero  at  all  solid 
and  fluid  boundaries.  The  boundary  conditions  are  updated  after  every  interior  points 
.updated  during  each  iteration.  Thus  the  boundary  values  as  well  as  interior  values  are 
iteratively  advanced  from  a  time  level  'n'  to  'n+T. 


Initial  Conditions  : 

In  the  case  of  external  flows,  we  assume  that  the  object  is  impulsively  started  from 
rest .  Therefore,  the  uniform  freestream  conditions  are  used  as  initial  conditions.  In  the  case 
of  internal  flows,  parallel  flow  solutions  (e.g.  Poiseulle  flow  in  a  square  duct)  are  used  to 
stan  the  calculations. 

Farfield  Boundary  Conditions  : 

For  external  flow  applications,  the  farfield  boundary  is  placed  far  away  from  the 
solid  surface.  Thus,  it  is  natural  to  specify  the  freestream  values  at  the  farfield  boundaries 
except  along  the  outflow  boundary  where  the  extrapolation  for  velocities  in  combination 
with  P  =  P «.  is  used,  to  account  for  the  removal  of  vorticity  from  the  flow  domain  by 
convective  process. 


Boundary  Conditions  on  the  Solid  Surface  : 

On  the  solid  surface,  the  no  slip  condition  is  imposed  for  velocity  components.  The 
surface  pressure  distribution  is  determined  by  solving  the  normal  gradient  of  pressure  to  be 
zero : 


(3.17) 


Some  researchers  (ref.22,  23)  obtain  the  boundary  conditions  for  pressure  from  the  normal 
component  of  momentum  equation  at  the  wall 


5p  _  1  d2un 
dn  Re  dn2 


(3.18) 


where  u  n  is  the  normal  component  of  velocity.  In  high  Reynolds  number  flows,  the 
viscous  stress  contribution  to  the  normal  momentum  equation  can  be  neglected  at  the  wall 
and  the  grid  point  adjacent  to  the  surface  will  be  sufficiently  fine  so  that  constant  pressure 
normal  to  the  surface  can  be  assumed.  Thus  Eq.(3.17)  is  an  acceptable  boundary  condition. 


Boundary  Conditions  on  the  Cut  and  Singular  Line  : 

Since  the  C-grid  and  the  H-0  grid  which  are  used  for  two-dimensional  airfoil 
problem  and  three-dimensional  body  of  revolution  have  a  cut  and  singularlines, 
respectively,  special  treatment  is  needed  (see  Fig.  3.2  and  3.3).  Across  the  cut  of  the  C- 
-grid  system,  flow  quantities  should  be  continuous.  The  flow  quantities  on  the  cut  can  be 
-obtained  by  averaging  the  flow  properties  from  above  and  below  the  cut.  On  the  singular 
lines  that  occur  in  a  H-O  grid  system,  the  flow  quantities  are  obtained  by  extrapolating  from 
two  adjacent  interior  points  and  then  averaging  them  azimuthally  to  ensure  that  the  flow 
quantities  are  singe-valued. 


Fig. 3. 3.  Singularline  of  the  H-O  grid  system 


3.5  Acceleration  by  Multigrid  Technique 

Since  the  matrix  [N]  (which  is  an  approximate  to  matrix  [M]  of  Eq.(3.12))  is  a 
simple  diagonal-matrix,  it  leads  to  slow  convergence  of  the  pressure  and  velocity  fields  at 
every  time  step.  Use  of  such  a  simple  diagonal  matrix  simplifies  the  inversion,  and  makes 
the  flow  solver  100%  vectorizable  and  parallelizable.  To  accelerate  the  present  procedure,  a 
multigrid  technique  (Coarse  Grid  Correction  method)  is  applied  in  this  study. 

The  principles  behind  the  present  multigrid  technique  are  as  follows.  The  quantities 
(Au,  Av,  Aw,  Ap)  may  be  viewed  as  Fourier  series-like  sums  made  of  components  of 
different  wave  lengths.  An  extremely  coarse  grid  linking  a  point  to  a  node  several  units 
away  is  effective  in  computing  the  long  wave  length  components.  A  very  fine  grid  is 
effective  in  computing  the  short  wave  length  components,  and  is  very  inefficient  for 
computing  the  long  wave  length  components.  The  multigrid  technique  attempts  to  compute 
these  individual  components  of  Aq  on  grids  of  several  levels  efficiendy.  When  the  process 
converges,  of  course,  the  discretized  equations  (i.e.  RHS  of  Eq.(3.7)  and  (3.11))  are 
exactly  satisfied  on  the  finest  grid. 

The  coarse  grid  correction  algorithm  presently  used  (given  here  for  2-grid  sequence 
for  simplicity)  is  as  follows  : 

i)  Compute  the  residual  {R}  appearing  on  the  right  hand  side  of  Eq.(3. 12)  on  the  fine  grid 

n+l.k 

using  q 

ii)  Transfer  the  residual  from  the  fine  grid  to  the  coarse  grid  using  the  injection  operation, 

T2hD 

1  h  K .  An  injection  operation  is  given  at  any  node  (i,j)  in  two-dimensional  case  by 


r  R  =  R  +  -rrf  R  .  +  R  .  .  +  R  +  R  .) 

h  ».J  1.J  2\  1+1-J  ‘-'-J  >.j+*  ‘.J-l/ 


(3.19) 


iii)  Compute  the  quantity  Aq  at  every  point  on  the  coarse  grid  by  solving  the  system  of 
equation  : 


[N]{Aq/J}  =  {lShR} 


(3.20) 


iv)  Interpolate  the  Aq  values  computed  in  step  (iii)  back  on  to  the  fine  grid  by  using  the 
bilinear  interpolation. 

v)  Compute  the  updated  values  of  the  flow  properties  qn+1,k+1  as  qn+1,k  +  Aq  . 


Repeat  step  (i)  -  (v)  till  Aq  is  driven  to  zero. 


The  present  2-D  solver  accepts  grids  upto  3  levels. 

To  the  writer's  knowledge,  the  multigrid  technique  in  unsteady  incompressible 
flows  has  been  applied  only  to  pressure-Poisson  equation.  The  u-,  v-  and  w-  momentum 
equations  are  usually  solved  only  on  a  single  grid.  The  present  work  fully  exploits  the 
-benefits  of  the  multigrid  method  for  all  the  equations,  while  keeping  the  form  of  the  matrix 
[N]  extremely  simple.  This  allows  use  of  larger  time  steps  and  improved  convergence  as 
discussed  on  Chapter  IV.  The  present  investigator  applied  a  conjugate  gradient  like  scheme, 
called  the  GMRES  (Generalized  Minimal  Residuals)  to  solve  Eq.(3.12).  The  matrix  [N] 
was  used  as  the  preconditioner.  The  success  of  the  GMRES  scheme  crucially  depends  on 

the  closeness  of  [N]  to  [M]  That  is  the  eigenvalues  of  the  matrix  [i  -N_1m|  must  be 

small  and  closely  packed.  The  use  of  GMRES  with  [N]  as  a  preconditioner  was  not 
successful. 


CHAPTER  IV 


RESULTS  AND  DISCUSSION 


In  this  chapter,  the  work  done  to  date  is  presented.  To  validate  the  present 
procedure,  three  cases  were  tested.  The  first  test  case  is  two-dimensional  unsteady  viscous 
flow  over  an  oscillating  airfoil.  The  second  is  three-dimensional  steady  flow  over  an 
ellipsoid  of  revolution.  The  third  is  the  flow  through  a  curved  duct.  Numerical  results  are 
presented  in  the  form  of  instantaneous  streamlines,  velocity  profiles,  vorticity  contours, 
surface  pressure  distribution,  and  aerodynamic  loads.  Streamlines  and  surface  pressure 
distributions  are  compared  with  flow  visualization  and  the  other  available  numerical  data  . 

4.1  Dynamic  Stall  of  an  Oscillating  Airfoil 

The  computations  were  carried  out  for  a  sinusoidally  pitching  NACA  0012  airfoil, 
at  Re  =  5,000  and  k  =  0.5,  where  k  is  reduced  frequency  of  oscillation,defined  by 


where  £2  is  the  radians  of  rotation  per  second  and  c  is  chord  of  airfoil.  The  physical 
interpretation  of  reduced  frequency  is  the  number  of  radians  of  oscillation  per  semi-chord 
length  of  travel.  This  case  has  been  previously  studied  by  Mehta  (ref.3)  at  NASA  Ames 
Research  Center  using  a  velocity-vorticity  formulation  and  its  fic>v  visualization  was 
carried  out  by  Werle  (ref.  24)  in  ONERA. 

After  the  flow  is  fully  developed  at  zero  angle  of  attack,  the  airfoil  is  allowed  to 
oscillate  in  pitch  through  an  angle  of  attack  range  from  0  degree  to  20  degree  given  by 
a  =  1 0°(  1  -  cos  t) .  Fig.l  shows  the  body-fitted  grid  around  the  airfoil  used  in  this  study. 
Fig. 2  shows  the  instantaneous  streamlines  (actually,  called  panicle  tracers  in  PLOT3D 
software),  velocity  profiles  and  vorticity  contours  at  selected  angle  of  attack.  Fig. 3  show  s 
the  surface  pressure  distribution.  In  general,  the  streamline  patterns  and  surface  pressure 
distributions  are  in  very  good  agreement  with  flow  visualization  and  Mehta's  numerical 


results  except  that  the  present  procedure  predicts  a  little  earlier  generation  of  vortex  than 
Mehta's  method.  The  flow  visualizations  were  carried  out  with  air  bubble  in  the  water 
tunnel.  Here,  we  should  note  that  photographs  showing  air  bubble  trajectories  were  taken 
at  an  exposure  time  of  1/10  seconds.  Therefore,  in  unsteady  flow  the  air  bubble  trajectories 
.near  the  surface  of  airfoil  represent  neither  streamlines  nor  streaklines  because  the  pictures 
contain  many  paths  over  the  exposure  time.  On  the  orther  hand,  the  instantaneous 
-streamline  is  a  streamline  at  any  instant  of  time,  i.e.  we  assume  the  flowfield  is  frozen  at 
any  instant  of  time  and  draw  the  streamline.  In  other  words,  the  instantaneous  streamline  is 
equivalent  to  the  bubble  trajectories  with  an  infinitesimal  exposure  time.  Thus,  the  flow 
visualization  with  air  bubble  is  different  from  the  instantaneous  streamline,  and  should  be 
used  only  for  qualitative  comparison.  Fig.4  shows  the  lift,  drag  and  moment  hysteresis 
loops.  The  main  feature  of  dynamic  stall  which  is  significantly  different  from  static  stall  is 
due  to  the  generation  of  a  vortex  near  the  leading  edge.  ThL  vor  »x  passes  over  the  upper 
surface  of  airfoil,  creating  large  variations  in  the  aerodynamic  forces  and  moment.  From 
these  figures,  it  is  seen  that  the  growth  of  lift  «i-‘ng  tl.e  upstroke  is  slow  and  gradual,  well 
past  the  static-stall  angle.  The  separation  region,  which  is  present  over  a  small  region  near 
the  trailing  edge  at  first,  moves  upstream  as  the  angle  of  attack  increases.  The  pitching 
moment  does  not  change  much  during  the  upstrok ,  The  surface  pressure  distribution  at  an 
angle  of  attack  of  18.6  degree  shown  in  Fig.3  shows  another  pressure  peak  near  the  quarter 
chord.  This  indicates  the  leading-edge  vortex  is  already  generated,  and  this  can  be  identified 
in  F*0.2  (c).  As  the  leading-edge  vortex  moves  downstream,  the  chordwise  surface 
pres  .’re  distribution  and  aerodynamic  forces  are  significantly  varied,  especially  during  the 
do'  -oke.  This  variation  may  depend  on  the  Reynolds  number,  airfoil  shapes  and 
reduced  frequency.  The  moment  stall,  associated  with  an  increase  of  negative  moment, 
begins  at  about  18.5  degree  in  the  downstroke. 


4.2  3-D  Steady  Flow  over  an  Ellipsoid  of  Revolution 

To  validate  the  capability  of  the  present  method  to  handle  three-dimensional  viscous 
flows,  the  present  procedure  was  tested  by  computing  the  flow  past  a  6:1  ellipsoid  of 
revolution  at  10  degree  angle  of  attack,  at  a  Reynolds  number  of  5,000.  Fig. 5  shows  the 
body-fitted  grid  system.  Fig. 6  shows  streamlines  over  the  body  surface.  There  is  a  limited 
amount  of  experimental  data  (ref.25,  26)  available  for  this  particular  configuration,  at  high 

A 

Reynolds  number  (Re=7.2  x  10  ).  Fig. 7  t  ;ws  the  surface  pressure  distribution  on  the 
windward  and  leeward  sides  of  the  symmetry  plane,  along  with  the  experimental  data. 


Good  agreement  is  evident  everywhere  except  in  the  last  10%  of  the  body,  where  the 
present  laminar  simulation  predicts  flow  separation,  and  a  flattening  out  of  the  pressure 
distribution. 


4.3  3-D  Steady  Flow  through  a  90°  Bended  Square  Duct 


To  validate  the  capability  to  handle  three-dimensional  internal  flow  problems,  the 
flow  within  a  square  duct  with  a  90-deg  bend  was  tested.  The  radius  of  curvature  of  the 
inner  wall  in  the  curved  section  is  1.8  times  of  the  side  length  of  square  cross-section.  This 
particular  configuration  (Fig. 8)  was  experimented  by  Humphrey  et  al.  (ref.27)  at  a 
Reynolds  number  of  790  based  on  the  average  inflow  velocity  and  hydraulic  diameter.  The 
inflow  and  outflow  velocity  profile  are  obtained  by  solving  the  equation  of  fully  developed 
duct  flow  (ref.28) : 


32u  82u 

dy*  +  a? 


1  dp 

=  — —  =  const. 
H  dx 


(4.2) 


This  equation  is  a  standard  form  of  Poisson  equation  and  can  be  solved  by  ADI  scheme. 
Fig.9  shows  the  streamwise  velocity  profiles  compared  with  the  experimental  data  of 
Humphrey  et  al.  The  plots  on  the  left  side  of  Fig.9  are  at  y/y  i /2  =  0-5,  which  is  midway 
between  the  left  side  wall  and  the  symmetry  plane  of  square  duct  and  the  right  side  plots  are 
at  y/y  1/2  =  0,  which  is  on  the  symmetry  plane.  The  inside  and  outside  wall  are 
corresponding  to  z  =  0  and  z  =  1,  respectively.  The  results  are  in  general  good  agreement 
with  experiments  except  at  6  =  90deg.  This  discrepancy  may  be  disappeared  if  a  more  fine 
grid  is  adopted.  This  will  be  further  investigated.  The  present  grid  system  is  61x21x21.  In 
Fig.  10,  the  cross-sectional  velocity  profiles  are  plotted  at  0  =  30,  60  and  90  deg.  The  top 
side  and  bottom  side  of  cross-section  are  the  inside  wall  and  outside  wall,  respectively.  In 
this  figure,  the  pair  of  secondary  vortices  is  shown  and  these  vortices  are  generated  due  to 
the  pressure  difference  between  the  higher  pressure  on  the  outside  wall  and  lower  pressure 
on  the  inside  wall.  Fig.l  1  is  a  sideview  of  streamwise  velocity  profiles  at  y/y1/2  =  0.5  and 
y/y  1/2  =  0-  Fig. 12  shows  streamwise  velocity  profiles  from  a  viewpoint  which  is  located  at 
upper  45°  in  the  xz-plane.  The  plot  at  z  =  0.25  is  corresponding  to  the  midway  plane 
between  the  inside  wall  and  the  plane  of  symmetry.  The  plots  at  z  =  0.5  and  0.75  are  on 
the  plane  of  symmetry  and  the  midway  between  the  outside  wall  and  the  symmetry  plane, 
respectively. 


4.4 


Acceleration  of  2-D  Flow  Solver  by  Multigrid  Technique 


The  multigrid  technique  was  implemented  to  the  two-dimensional  steady  and 
unsteady  solver.  The  fine  grid  system  has  (81x41)  grid  points  and  the  coarse  grid  system 
*has  the  half  of  the  fine  grid  points,  i.e.  (41x21)  grid  points,  and  the  coarsest  grid  system 
has  (21x11)  grid  points.  The  two  grid  system  consists  of  the  fine  and  coarse  grid  system 
(Fig.  4.1. (a))  and  the  three  grid  system  consists  of  all  of  them  as  shown  in  Fig.4.1.(b). 
Especially,  three  grid  system  such  as  Fig.  4.1.(b)  is  called  V-cycle. 


Coarsest  Grid 

(a) 


(b) 


Fig.  4.1  Structure  of  multigrid  cycle 

(a)  Two-grid  system 

(b)  Three-grid  system 


Fig.  13  shows  the  convergence  history  of  the  global  residual  (l2-norm  of  RHS  of  Eq 

(3.12))  reduction  in  CPU  time  for  steady  flow  over  NACA  0012  airfoil  at  zero  angle  of 
attack.  Upto  40%  and  60%  acceleration  was  obtained  using  two-  and  three-grid  system, 
respectively.  The  CPU  time  is  based  on  25  iterations  at  each  time  step  on  an  IBM 
RISC/6000  workstation.  Fig.  14  shows  the  history  of  global  residual  for  sinusoidally 
oscillating  airfoil  (50  iterations/time  step),  where  the  three-grid  system  is  used  for 
multigrid.  The  residual  by  the  multigrid  technique  maintains  lower  level  than  that  of  single 
grid  iteration  procedure  indicating  that  the  discretized  equations  are  solved  to  much  high 
levels  of  accuracy  using  the  multigrid  technique.  The  surface  pressure  distribution  and 


dynamic  stall  hysteresis  is  nearly  the  same  as  those  of  single  grid  system  and  are  not  plotted 
here. 
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CHAPTER  V 

PROPOSED  WORK 


An  iterative  procedure  for  two-  and  three-dimensional  unsteady,  incompressible,  viscous 
flow  has  been  developed.  It  has  been  applied  to  massively  separated  flow  over  oscillating  airfoil, 
three-dimensional  flow  past  an  ellipsoid  of  revolution,  and  three  dimensional  flow  through  a 
curved  square  duct.  Good  agreement  with  published  experimental  and  numerical  data  has  been 
obtained.  After  the  validation  of  the  present  procedure,  techniques  for  acceleration  were  explored. 
It  was  found  that  the  multigrid  technique  was  efficient  in  reducing  the  CPU  time  needed  for  the 
simulation  and  improved  the  solution  quality  because  of  the  lower  residuals  achieved.  The  GMRES 
does  not  work  successfully  presumably  because  of  the  diagonal  algorithm  used  as  the 
preconditioner. 

The  present  multigrid  iterative  scheme  for  unsteady  incompressible  viscous  flows  is  being 
extended  to  three-dimensions.  The  single-grid  version  has  already  been  tested  as  discussed  in 
Chapter  IV.  It  is  proposed  that  the  following  calculations  be  done  to  test  the  suitability  of  this 
procedure  to  three-dimensional  viscous  flows  : 

(a)  Completion  of  three-dimensional  multigrid  solver  validation  for  a  curved  square  duct  flow 

(b)  Application  of  the  three-dimensional  multi  grid  method  to  flow  past  a  marine  propeller  (Fig.  15) 
or  a  direct  numerical  simulation  of  the  turbulent  Rayleigh-Benard  problem 

(c)  Extension  to  higher  order  accuracy  in  time.  Specifically,  the  following  studies  will  be  done  : 
This  is  simple,  and  simply  requires  replacement  of  terms  such  as  (  qn+1' k  —  qn  )  /  At  appearing 
in  Eq.(3.9)  to  (  3  qn+1,k  -  4  qn  +  qn_1  )  /  At  etc.  One  can  show  this  formally  raises  the  time 
accuracy  to  (At)2. 


(b) 


Figure  1.  Body-Fitted  Grid  Around  a  NACA  0012  airfoil 


(a)  a  =  0  (deg) 

Figure  2.  Instantaneous  Streamlines,  Velocity  Profiles,  and  Vorticity 
Contours  at  Selected  Angle  of  Attack  with  Experimental  Flow 
Visualizations 


(b)  a  =  14.6  (deg) 


Figure  2.  Continued.. 


(d)  a  =  20  (deg) 
Figure  2.  Continued. 


(e)  a  =  19.8  (deg) 
Figure  2.  Continued. 


Figure  2.  Continued. 
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Figure  3.  Continued 
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Figure  5.  Body-Fitted  Grid  Around  an  Ellipsoid  of  Revolution 


Figure  8.  Body-fitted  grid  within  a  square  duct 
with  a  90-deg  bend. 


Figure  9.  Streamwise  velocity  profiles  compared  with  experiments. 


(b)  0  =  60 


/  _ 


_  \ 


\ 


\ 


"  \  \  *//" 

A  1  i  '  /  /  ' 


! 


'll1  \  \  ' 
/  /  1  '  1  \  \ 
/  1  1  *  l  \  \ 
//Ui|\ 


/ 


// 1 


(c)  0  =  90 


Figure  10.  Cross-sectional  velocity  profiles  at  three 
streamwise  stations  in  the  curved  section. 


Figure  11.  Streamwise  velocity  profiles  at  y/y1/2=  0.5  and  y/y1/2=  0- 
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Figure  13.  Comparison  of  Convergence  History  of  NACA  0012  Airfiol  for 
Steady  State  ( at  Zero  Degree  Angle  of  Attack  ) 
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Figure  14.  Comparison  of  Global  Residual  History  of  a  Sinusoidally 
Oscillating  Airfoil 
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